Otimização

Otimização#

Considere o seguinte problema. Desejamos construir uma calha a partir de uma chapa de metal com \(30\) cm de largura em forma de trapézio conforme a figura abaixo.

Hide code cell source
# HIDE CODE
import plotly.express as px
import plotly.graph_objects as go
import numpy as np
import sympy as sp
Hide code cell source
# HIDE CODE
theta = np.linspace(0, np.pi/3,30)

fig = go.Figure(
    data=[go.Scatter(x=[-0.5 + np.cos(2*np.pi/3), -0.5, 0.5, 0.5 + np.cos(np.pi/3)], 
                     y=[np.sin(2*np.pi/3), 0.0 , 0.0 , np.sin(np.pi/3)],
                     mode="lines",
                     showlegend=False,
                     line=dict(width=4, color="blue")),
        go.Scatter(
    x=[-0.45 + np.cos(2*np.pi/3)/2, 0, 0.45 + np.cos(np.pi/3)/2,0.625],
    y=[np.sin(2*np.pi/3)/2 + 0.05, -0.05, np.sin(np.pi/3)/2+ 0.05,0.15],
    mode="text",
    text=[r"$x$", r"$30-2x$", r'$x$', r'$\theta$'],
    textposition="bottom center",
    showlegend=False,
    textfont=dict(
        #family="sans serif",
        size=20,
        color="black"
    )
                ),
        go.Scatter(x=0.5 + 0.1*np.cos(theta), 
                     y=0.1*np.sin(theta),
                     mode="lines",
                     showlegend=False,
                     line=dict(width=2, color="black")),
    ]
)

fig.add_shape(type="line",
    x0=0.5, y0=0, x1=1, y1=0,
    line=dict(
        color="blue",
        width=4,
        dash="dot",
    )
)

fig.update_yaxes(
    scaleanchor = "x",
    scaleratio = 1,
  )

fig.update_layout(xaxis=dict(range=[-1,1],showticklabels=False),yaxis=dict(range=[-0.3,1],showticklabels=False))



fig.show()

Função área da seção transversal

\(A(x, \theta) = (2(30 - 2x)+ 2x \cos(\theta) )\dfrac{x \sin(\theta)}{2} \)

ou seja

\(A(x, \theta) = (30 - 2x+ x \cos(\theta) )x \sin(\theta) \)

domínio

\( \mathcal{D} = \left\\{ (x,\theta) \in \mathbb{R}^2 \ \ | \ \ 0 \leq x \leq 15 \quad \text{e} \quad 0 \leq \theta \leq \dfrac{\pi}{2} \right\\} \)

Derivadas parciais

\(\dfrac{\partial A(x,\theta)}{\partial x} = (\cos(\theta) - 2)x\sin(\theta) + (30 - 2x+ x \cos(\theta) )\sin(\theta)\)

\(\dfrac{\partial A(x,\theta)}{\partial \theta} = -x^2\sin^2(\theta) + (30 - 2x+ x \cos(\theta) )x\cos(\theta)\)

Queremos encontrar os pontos críticos no interior do domínio

\[\begin{equation*} \begin{cases} \dfrac{\partial A(x,\theta)}{\partial x} = (\cos(\theta) - 2)x\sin(\theta) + (30 - 2x+ x \cos(\theta) )\sin(\theta) &= 0 \\ \dfrac{\partial A(x,\theta)}{\partial \theta} = -x^2\sin^2(\theta) + (30 - 2x+ x \cos(\theta) )x\cos(\theta) &= 0 \end{cases} \end{equation*}\]

Resolvemos o sistema acima, multiplicando a primeira equação por \(x \cos(\theta)\) e subtraimos a segunda equação multiplicada por \(\sin(\theta)\), de onde obtemos

\[\begin{align*} x^2 \sin(\theta) \left( (\cos(\theta) - 2)\cos(\theta) + \sin^2(\theta)\right) &= 0 \quad 0 < x < 30, \quad 0 < \theta < \frac{\pi}{2}\\ \cos^2(\theta) - 2\cos(\theta) + 1 - \cos^2(\theta)&= 0 \\ - 2\cos(\theta) + 1 &= 0 \\ \cos(\theta) &= \frac{1}{2} \quad 0 < \theta < \frac{\pi}{2} \\ \theta & = \frac{\pi}{3} \end{align*}\]

Substituindo o valor de \(\theta = \frac{\pi}{3}\) em qualquer uma das equações do sistema acima, encontramos \(x=10\).

Hide code cell source
# HIDE CODE
points = 100

x = np.linspace(0, 15, points)
theta = np.linspace(0, np.pi/2, points)

xGrid, thetaGrid = np.meshgrid(x, theta)

z = (30 - 2*xGrid + xGrid*np.cos(thetaGrid))*xGrid*np.sin(thetaGrid)

fig = go.Figure(data=[go.Surface(x=x, y=theta,z=z)])

fig.show()

Observe que nesse caso o valor máximo absoluto é atingido no interior do domínio e não da fronteira.

As curvas de nível#

Hide code cell source
# HIDE CODE
points = 100

x = np.linspace(0, 15, points)
theta = np.linspace(0, np.pi/2, points)

xGrid, thetaGrid = np.meshgrid(x, theta)

z = (30 - 2*xGrid + xGrid*np.cos(thetaGrid))*xGrid*np.sin(thetaGrid)

fig = go.Figure(data=[go.Contour(x=x, y=theta,z=z)])

fig.show()
Hide code cell source
# HIDE CODE
points = 100

x = np.linspace(0, 15, points)
theta = np.linspace(0, np.pi/2, points)

xGrid, thetaGrid = np.meshgrid(x, theta)

z = (30 - 2*xGrid + xGrid*np.cos(thetaGrid))*xGrid*np.sin(thetaGrid)

fig = go.Figure(data=[go.Surface(x=x, y=theta,z=z)])

fig.update_traces(contours_z=dict(show=True, usecolormap=True,
                                  highlightcolor="limegreen", project_z=True))

fig.show()

Cálculo com sympy#

Primeiro vamos usar as variáveis \(x, \theta\) como símbolos para que o sympy funcione corretamente.

x, theta = sp.symbols('x theta')

Vamos definir a função \(A(x, \theta)\)

def A(x,theta):
  return (30 - 2*x + x*sp.cos(theta))*x*sp.sin(theta)

Calcular as derivadas de \(A\), primeiro em relação a \(x\).

A_x = sp.diff(A(x,theta),x)
A_x
\[\displaystyle x \left(\cos{\left(\theta \right)} - 2\right) \sin{\left(\theta \right)} + \left(x \cos{\left(\theta \right)} - 2 x + 30\right) \sin{\left(\theta \right)}\]

Agora em relação a \(\theta\)

A_theta = sp.diff(A(x,theta),theta)
A_theta
\[\displaystyle - x^{2} \sin^{2}{\left(\theta \right)} + x \left(x \cos{\left(\theta \right)} - 2 x + 30\right) \cos{\left(\theta \right)}\]

Perceba que resolvemos o exercício acima com a derivada calculada corretamente segundo o sympy.

Vejamos como utilizar o sympy agora para encontrar os pontos críticos.

sp.solve([A_x, A_theta], [x,theta], dict=True)
[{theta: 0, x: 0},
 {theta: 0, x: 30},
 {theta: -pi/3, x: 10},
 {theta: pi/3, x: 10}]

Observe que apenas o ponto \((x,\theta) = (10,\frac{\pi}{3})\) pertence ao interior do domínio \(\mathcal{D}\).

O exercício \(3\) da seção “Problemas Quentes” do capítulo \(14\) de [SCW22] pergunta se seria melhor dobrar a folha de metal para fazer uma calha com seção tranversal semicircular. Vamos comparar os dois modelos.

O valor da área da seção transversal em formato de trapézio dado pela solução que encontramos anteriormente é:

Utilizando sympy.

A(10,sp.pi/3)
\[\displaystyle 75 \sqrt{3}\]

Utilizando o numpy.

A(10,np.pi/3)
\[\displaystyle 129.903810567666\]

Uma outra maneira de verificar a expansão decimal usando o sympy seria

sp.N(A(10,sp.pi/3))
\[\displaystyle 129.903810567666\]

Agora vejamos o valor da área com seção tranversal semicircular, primeiro vamos calcular o raio.

\[\begin{equation*} \pi r = 30 \Longrightarrow r = \frac{30}{\pi} \end{equation*}\]

Portanto a área da seção tranversal semicircular é

\[\begin{align*} \frac{1}{2} \pi r^2 &= \frac{1}{2} \pi \left( \frac{30}{\pi} \right)^2 \\ &= \frac{450}{\pi} \end{align*}\]
450/np.pi
143.2394487827058

Ou seja, a calha semicircular tem uma capacidade maior.

Vejamos um gráfico comparando as duas figuras.

Hide code cell source
# HIDE CODE
raio = 3/np.pi
theta = np.linspace(-np.pi, np.pi,300)

fig = go.Figure(
    data=[go.Scatter(x=[-0.5 + np.cos(2*np.pi/3), -0.5, 0.5, 0.5 + np.cos(np.pi/3)], 
                     y=[np.sin(2*np.pi/3), 0.0 , 0.0 , np.sin(np.pi/3)],
                     mode="lines",
                     showlegend=False,
                     line=dict(width=4, color="blue")),
        go.Scatter(x=raio*np.cos(theta), 
                     y=raio + raio*np.sin(theta),
                     mode="lines",
                     showlegend=False,
                     line=dict(width=4, color="red")),
    ]
)

fig.update_yaxes(
    scaleanchor = "x",
    scaleratio = 1,
  )

fig.update_layout(xaxis=dict(range=[-1,1],showticklabels=False),yaxis=dict(range=[-0.3,1],showticklabels=False))



fig.show()